EcoCommons Marxan MaPP connection

EcoCommons -> Marxan MaPP connection

Author: “Author Name”

Date: “2024-09-10”

Introduction

Using the Species distribution modeling techniques provided by the EcoCommons Platform (www.ecocommons.org.au), we produced probability distribution maps for the three Queensland endangered species: koala, brush tailed rock-wallaby, and beach stone curlew.

Then we adjusted the probability distribution maps of these three species with the planning units shapefile prepared by the Marxan MaPP, and ran four planning scenarios with a target of expanding the coverage of protected areas in QLD to 30%.

EcoCommons Outputs

Species records pulled from GBIF, ALA, EcoPlots, OBIS
Species distribution modelling output: Species distribution Probability maps (This is the input tested in this project).

Marxan MaPP Inputs

Shapefile of planning area and units.
Shapefile of cost.
Shapefile and csv of biodiversity features (Where EcoCommons can help!).

EcoCommons connects with Marxan Showcase:

Make sure you are in the directory you want

getwd()

[1] "/Users/zhaoxiang/Documents/tmp/ecocommons-marxan/ecocommons-marxan-integration-poc"

Install and load essential packages

# First specify the packages of interest
packages = c("sf", "terra", "ggplot2", "ggspatial", "raster", "dplyr")

# Now load or install&load all. This process will take a long time since we are using a virtual environment and install a lot of packages.
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
     }
  }
 )

1. We get the QLD planning units from Marxan MaPP

QLD_Unit <- "qld_3species_Marxan/QLD_plannningunits/cost-surface-template_dbfd87df-99bc-4229-a607-04e632a67fda_20240828T015655.shp"  #This cost-surface-template was prepared by the Marxan Mapp with a resolution of 189 Km2, which is the highest resolution Marxan Mapp can give at this scale.

QLD_Unit  <- st_read(QLD_Unit)

# Calculate the resolution since Marxan MaPP for visulization purpose
areas <- st_area(QLD_Unit)
areas_numeric <- as.numeric(areas)
average_area <- mean(areas_numeric)
# Convert to numeric
average_area_km2 <- average_area / 1e6

# Get the number of rows
n_rows <- nrow(QLD_Unit)

# Plot the shapefile with no fill color and number of rows in the title
ggplot(data = QLD_Unit) +
  geom_sf(fill = NA, color = "gray") +
  theme_minimal() +
  ggtitle(paste("QLD Planning Units:", n_rows, "\n",
                "Resolution of planning in square kilometers:", round(average_area_km2)))+
  theme(plot.title = element_text(hjust = 0.5))  # Center the title

2. I made a cost layer using the reciprocal of the distance to state-owned road as a surrogate of the cost.

The assumption is: the closer to the state owned road, the more expensive to purchase the unit.

QLD_cost_road <- st_read("qld_3species_Marxan/QLD_Cost/QLD_cost_road.shp")

# Plot the shapefile with continuous cost_road values
ggplot(QLD_cost_road) +
  geom_sf(aes(fill = cost_road)) +
  scale_fill_continuous(name = "Cost",
                        low = "lightblue", high = "red",
                        labels = c("0 (Low cost)", "1 (High cost)"),
                        breaks = c(0.01, 1)) +
  theme_minimal() +
  labs(title = "Cost: using the distance to road of each Unit as a proxy")+
  theme(plot.title = element_text(hjust = 0.5))  # Center the title

3. Biodiversity features. I used EcoCommons to produce three species’ SDM to start with.

Species 1: koala
Species 2: brush tailed rock-wallaby
Species 3: beach stone curlew

# Define the folder path where the rasters are stored
folder_path <- "qld_3species_Marxan/QLD_feature/"

# Get a list of all .tif files in the folder
raster_files <- list.files(path = folder_path, pattern = "\\.tif$", full.names = TRUE)

# Extract the species names from the file names (removing the folder path and .tif extension)
species_names <- tools::file_path_sans_ext(basename(raster_files))


# Read all raster files in one go using lapply
raster_list <- lapply(raster_files, rast)  # Use rast() from terra for reading rasters

# Using QLD_Unit as the spatial vector for masking

# Transform the raster CRS to match the vector CRS and apply masking in one step
raster_list <- lapply(raster_list, function(r) {
  r_transformed <- project(r, crs(vect(QLD_Unit)))
  mask(r_transformed, vect(QLD_Unit))
})

# Function to convert rasters to data frames and combine them
prepare_raster_data <- function(raster_list, species_names) {

  # Initialize an empty data frame
  combined_df <- data.frame()
  # Loop through each raster and combine them into one data frame
  for (i in seq_along(raster_list)) {
    # Convert raster to a data frame
    raster_df <- as.data.frame(raster_list[[i]], xy = TRUE)
    # Rename the third column to 'value' or any appropriate name for the raster values
    names(raster_df)[3] <- "value"
    # Add a column to identify the species name
    raster_df$species <- species_names[i]
    # Combine the raster data with the overall data frame
    combined_df <- bind_rows(combined_df, raster_df)
}
  return(combined_df)
}

# Prepare the combined data frame
combined_raster_df <- prepare_raster_data(raster_list, species_names)

# Create the ggplot with facet_wrap to display each raster in a separate facet
ggplot(combined_raster_df, aes(x = x, y = y, fill = value)) +  # Use the correct column name for fill
  geom_raster()+
  facet_wrap(~ species, ncol = 3) +  # Adjust ncol to control the number of columns
  scale_fill_viridis_c() +  # You can adjust the color scale as needed
  labs(title = "Species SDM") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

4. We need to turn these SDMs to binary results (shapefies).

 # Define a function to process each species
 process_species <- function(raster_data, QLD_Unit, species_name, output_dir) {
  # Check the input raster and vector objects
  print(raster_data)  # Check raster structure
  print(QLD_Unit)     # Check vector structure

  # Reproject raster to match vector CRS
  raster_data_transformed <- project(raster_data, crs(vect(QLD_Unit)))
  print(raster_data_transformed)  # Check the transformed raster

  # Extract mean raster values for each polygon
  extracted_values <- extract(raster_data_transformed, vect(QLD_Unit), fun = mean, na.rm = TRUE)
  print(extracted_values)  # Debugging: check extracted values

  # Update the QLD_Unit feature column
  names(QLD_Unit)[names(QLD_Unit) == "cost"] <- "feature"
  QLD_Unit$feature <- extracted_values[, 2]

  # Subset the polygons where the feature value is >= 0.5
  QLD_species <- subset(QLD_Unit, feature >= 0.5)
  print(QLD_species)  # Check the subset for debugging

  # Write the subset to a shapefile
  shapefile_base <- file.path(output_dir, species_name)
  st_write(QLD_species, paste0(shapefile_base, ".shp"), delete_layer = TRUE)

  # Zip the shapefile components
  shapefile_files <- list.files(output_dir, pattern = paste0(species_name, "\\."), full.names = TRUE)
  shapefile_base_names <- basename(shapefile_files)
  zipfile_path <- file.path(output_dir, paste0(species_name, ".zip"))
  old_wd <- setwd(output_dir)
  zip(zipfile_path, shapefile_base_names)
  setwd(old_wd)
}

# Set the output directory
output_dir <- "/Users/zhaoxiang/Documents/tmp/ecocommons-marxan/ecocommons-marxan-integration-poc/qld_3species_Marxan/QLD_feature/Marxan_feature_input"  #change these to the absolute path to the folder you would like to save your Marxan MaPP Feature input

# Apply the function to each species in raster_list
for (i in seq_along(raster_list)) {
  process_species(raster_list[[i]], QLD_Unit, species_names[i], output_dir)
}

5. Plot species SDM binary shapefile outputs for double check

# List all the shapefiles in the directory (assuming each species has its own shapefile)
species_files <- list.files(output_dir, pattern = "\\.shp$", full.names = TRUE)

# Extract species names from the filenames (you can adjust this depending on your naming conventions)
species_names <- tools::file_path_sans_ext(basename(species_files))

# Load all species shapefiles and add a species identifier
species_sf_list <- lapply(seq_along(species_files), function(i) {
  sf <- st_read(species_files[i])
  sf$species <- species_names[i]  # Add species name column
  return(sf)
})

# Combine all species into one dataset
combined_species_sf <- do.call(rbind, species_sf_list)

# Plot the unit (base map) first and overlay the species habitats without borders
combined_plot_with_unit <- ggplot() +
  geom_sf(data = QLD_Unit, fill = NA, color = "grey", size = 0.01) +  # Base map (QLD Unit)
  geom_sf(data = combined_species_sf, aes(fill = species), color = NA) +  # No borders for species
  scale_fill_manual(values = RColorBrewer::brewer.pal(n = length(species_names), name = "Set1")) +  # Automatically assign colors
  theme_minimal() +
  labs(title = "Species Habitats within QLD Unit",
       subtitle = paste(species_names, collapse = ", ")) +  # List all species in subtitle
  theme(legend.title = element_blank())

# Display the plot
print(combined_plot_with_unit)

6. Plot species SDM binary shapefile outputs for double check

# Function to extract presence (1) and absence (0) from raster based on a threshold (e.g., 0.5)

extract_presence_absence <- function(raster_data, unit) {
  extracted_values <- extract(raster_data, vect(unit), fun = mean, na.rm = TRUE)
  presence_absence <- ifelse(extracted_values[, 2] >= 0.5, 1, 0)
  return(presence_absence)
}

# Create an empty presence-absence data frame
presence_absence_df <- data.frame(puid = QLD_Unit$puid)  # Assuming 'puid' is the unique identifier

# Loop through each species raster in the raster list and extract presence-absence data
for (i in seq_along(raster_list)) {
  # Generate a dynamic presence column name for the current species
  presence_col_name <- paste0(species_names[i], "_presence")
  
  # Extract presence/absence data and add it to the presence-absence dataframe
  presence_absence_df[[species_names[i]]] <- extract_presence_absence(raster_list[[i]], QLD_Unit)
}

# Write the final presence-absence data frame to a CSV file
output_csv <- file.path(output_dir, "presence_absence_species.csv")
write.csv(presence_absence_df, output_csv, row.names = FALSE)

# Check the CSV output
print(head(presence_absence_df))

  puid beach_stone_curlew_GLM brushtailed_rockwallaby_GLM Koala_GLM
1    1                      0                           0         0
2    2                      0                           0         0
3    3                      0                           0         0
4    4                      0                           0         0
5    5                      0                           0         0
6    6                      0                           0         0

Marxan Four scenarios solutions:

Our SDMs input to Marxan MaPP:

EcoCommons SDMs output of three species on Marxan MaPP

Scenario 1: No SDMs, No Costs